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PREFACE 


This  is  a  supplement  to  the  Harwell  Subroutine  Library's  catalogue  R.7477  and  gives  a  list  of 
the  subroutines  introduced  into  the  library  in  the  period  from  August  1974  to  July  1977.  We 
have  not  included  a  general  index  to  the  the  library  this  time.  The  total  list  for  the  whole  library 
is  now  covered  by  the  report  R.7477  and  its  two  supplements  (1974)  and  (1977)  plus  the  deletion 
of  the  subroutines  EB01A/AD,  EB04A/AD,  LAO  1  A/AD,  MCI  3  A/C,  OEOIA/B,  PA01A/AD, 
PA02A/AD,  and  TG01A/AD. 

Many  of  our  new  subroutines  have  been  written  with  portability  in  mind  and  it  is  the  intention 
of  some  of  our  authors  to  write  code  that  conforms  to  the  Fortran  IV  ANSI  standard  (1966).  To 
do  this  we  use  the  Bell  Laboratories’  verifier  program,  described  in  Ryder,  ‘The  PFORT  Verifier’, 
Software  Practice  and  Experience,  4,  (1974),  to  check  out  the  new  subroutines.  We  have 
indicated  those  subroutines  which  have  passed  the  verifier  test  by  specifying  the  language  as 
‘Fortan  IV  (standard)’.  There  are  many  other  new  subroutines  that  we  are  not  flagged  as  verified 
which  are  virtually  standard  but  have  failed  the  test  on  some  small  point,  e.g.  OE02A  is  standard 
except  for  a  single  use  of  the  IBM  END=  parameter  in  a  READ  statement,  and  in  many  cases  a 
standard  version  can  be  obtained  using  OE04A  to  process  the  source  deck  and  these  we  have 
marked  ‘(standard  available)’. 

Again,  with  portability  in  mind,  we  have  begun  to  include  Fortran  versions  of  our  machine 
language  subroutines  on  the  source  tape  and  these  we  have  flagged  by  specifying  the  language  as 
370/BAL  (Fortran  available). /\ 

Harwell  Subroutine  Librarian 
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EA12A 


FM02AS 


Given  a  real  symmetric  nutria  A  =  {a^ }  of  order  n, 
calculates  the  k  largest  eigenvalues,  A  j>A2  A^ 

and  their  associated  eigenvectors  x(  /=  1,2,... Jc,  where 

A*/ -Vi  • 

The  subroutine  uses  the  method  of  simultaneous 
iteration  and  also  allows  the  user  to  take  advantage  of 
sparsity.  Eigensoiutions  from  other  parts  of  the 
spectrum  can  be  obtained  by  using  shifts  of  the  form 
A  —  0 1  optionally  combined  with  inversion. 

Versions:  EA12A,  EA12AD. 

Calls:  EA13A,  FA01AS,  FM01AS,  FM02AS. 
Language:  Fortran  IV  (standard). 

Date:  June  1976. 

Size:  8. IK;  932  cards. 

Origin:  I.S.DufT,  A.E.R.E.,  Harwell. 


EA13A 

Given  a  real  symmetric  matrix  A  ,  dads  al  its 
eigenvalues  A,  and  eigenvectors  x/ l=l,2,...,n,  where 

Axt  =  A^  . 

The  classical  Jacobi  method  is  used. 

Remark:  Supersedes  EA03A. 

Versions:  EA13A,  EA13AD. 

Language:  Fortran  IV  (standard). 

Date:  June  1976. 

Size:  3.3K;  168  cards. 

Origin:  I.S.DufT,  A.E.R.E.,  Harwell. 


FM01AS 

To  evaluate  the  inner  product  of  two  real  vectors, 

given  vectors  a »  {a/ }  and  b  =  {i>/ }  of  length  «, 
evaluates 


To  evaluate  the  inner  product  of  two  real  vectors 
when  the  elements  of  each  vector  are  stored  at  some 
fixed  displacement  from  neighbouring  elements.  Given 
vectors  a  =  {a/ }  and  b  =  {bt }  of  length  n  evaluates 

n 

W=1  aA 

/=i 

The  subroutine  can  be  used  to  evaluate  inner 
products  involving  rows  of  multi-dimensioned  arrays. 
It  is  a  FUNCTION  subroutine  written  in  machine 
language  and  is  fast  in  execution. 

Remark:  Supersedes  MC03AS. 

Versions:  FM02AS,  FM02AD,  FM02AQ.* 
language:  370/BAL  (Fortran  available). 

Date:  April  1975. 

Size:. IK;  125  cards. 

Origin:  MJ.Hopper,  A.E.R.E.,  Harwell. 


FM03AS 

To  evaluate  the  triple  inner  product  of  three  real 
vectors,  given  vectors  a  =  ia()  ,  b  =  {6,}  and  c  =  {c(( 
of  length  n,  evaluates 

n 

w=y  o.b.c. 

1=1  1 

The  subroutine  is  a  FUNCTION  routine  written  in 
machine  code  and  is  fast  in  execution. 

Remark:  Supersedes  MC05AS. 

Versions:  FM03AS,  FM03AD,  FM03AQ* 
Language:  370/BAL  (Fortran  available). 

Date:  May  1975. 

Size:  .IK;  115  cards. 

Origin:  MJ.Hopper,  A.E.R.E.,  Harwell. 


n 


FM04AS 


The  subroutine  is  a  FUNCTION  routine  written  in 
machine  code  and  it  fast  in  execution. 

Remark:  Supersedes  MC02AS. 

Versions:  FM01AS,  FM01AD,  FM01AQ  * 
Language:  370/BAL  (Fortran  available). 

Date:  January  1975. 

Size:  .IK;  105  cards. 

Origin:  MJ.Hopper,  A.E.R.E.,  Harwell. 


To  evaluate  the  triple  Inner  product  of  three  real 
vectors  when  the  elements  of  each  vector  are  stored  at 
some  fixed  displacement  from  neighbouring  elements. 
Given  vectors  a  =  (a(),b  =  (i(|  and  c  =  {c,}  of  length 
n,  evaluates 


w 


-  S«,  V, 


Q  versions  give  extended  precision. 


The  subroutine  is  a  FUNCTION  routine  written  in 
machine  code  and  is  fast  in  execution. 


It  is  a  FUNCTION  subroutine  written  in  machine 
language  and  is  fast  in  execution. 


Versions:  FM04AS,  FM04AD,  FM04AQ* 
Language:  370/BAL  (Fortran  available). 
Date:  May  1975. 

Sise:  .2K;  148  cards. 

Origin:  MJ.Hopper,  A.E.R.E.,  Harwell. 


Remark:  Supersedes  ME06AS. 

Versions:  FM06AS,  FM06AD,  FM06AQ* 
Language:  370/BAL  (Fortran  available). 
Date:  June  1975. 

Size:  .2K;  180  cards. 

Origin:  MJ.Hopper,  A.E.R.E.,  Harwell. 


FMOSAS 

To  evaluate  the  inner  product  of  two  complex 
vectors,  given  complex  vectors  a  =  (  and  b  =  {6( }  of 

length  n,  evaluates 

it 

i=\ 

or  optionally  evaluates 


it 


where  the  b(s  are  the  complex  conjugates  of  the  i>(s. 

The  subroutine  is  a  FUNCTION  routine  written  in 
machine  code  and  is  fast  in  execution. 

Versions:  FM05AS,  FM05AD,  FM05AQ.* 
Language:  370/BAL  (Fortran  available). 

Date:  June  1975. 

Size:  -2K;  150  cards. 

Origin:  MJ.Hopper,  A.E.R.E.,  Harwell. 


FM06AS 

To  evaluate  the  inner  product  of  two  complex 
vectors  when  the  elements  of  each  vector  are  stored  at 
some  fixed  displacement  from  neighbouring  elements. 
Given  complex  vectors  a  =  ia,}  and  b  =  {bt }  of  length 
n,  evaluates 


n 


or  optionally  evaluates 

n 


where  the  6f  s  are  the  complex  conjugates  of  the  b{ s. 

The  subroutine  can  be  used  to  evaluate  inner 
products  involving  rows  of  multi-dimensioned  arrays. 


*  Q  versions  give  extended  precision. 


LA01B 

Solves  the  linear  programming  problem,  finds 
x  =  \xj  )n  which  minimizes  the  linear  function 

it 

/to  =  2  c,x, 

/=!  J  J 

subject  to  the  linear  constraints 

n 

2  ayxj  ^  bt  1=1,2,... yk 

J=\  v  1 

It 

£a.,x.=b.  l=fc+l . m 

J=  \  u  J 

where  1=1,2,. ...n. 

The  Revised  Simplex  method  is  used  where  an 
inverse  of  the  basis  is  maintained  and  updated  at  each 
iteration.  There  is  an  option  to  allow  the  basis  inverse 
to  be  cleaned  up  at  the  end  of  the  calculation  and  the 
solution  checked  for  ill-conditioning  and  if  necessary 
the  calculation  is  restarted.  Another  option  allows  the 
user  to  force  specific  columns  of  the  LP  tableau  into 
the  initial  basis. 

LAO  IB  returns  an  error  indicator  and  an  estimate  of 
the  ill-conditioning.  Several  print  options  are  offered. 

Remark:  Supersedes  LA01A. 

Versions:  LA01B,  LA01BD. 

Calls:  FM01AS. 

Language:  Fortran  IV. 

Date:  January  1975. 

Size:  8.1K;  684  cards. 

Origin:  MJ.Hopper,  A.E.R.E.,  Harwell. 


LA05A 

To  factorize  a  matrix,  solve  corresponding  systems 
of  linear  equations  and  update  the  factorization  when  a 
column  of  the  matrix  is  altered,  exploiting  sparsity  in 
all  cases.  Its  primary  application  is  likely  to  be  for 
handling  linear  programming  bases. 
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The  subroutine  decomposes  the  matrix  into 
triangular  factors  using  sparse  matrix  techniques 
similar  to  those  described  in  Curtis  and  Reid,  Harwell 
report  R.6844,  (1971). 

Versions:  LAOSA,  LAOSAD. 

Calls:  MC20A. 

Language:  Fortran  IV. 

Date:  December  197S. 

Size:  20.6K;  959  cards. 

Origin:  J.K.Reid,  A.E.R.E.,  Harwell. 


hand  sides  to  solve  further  equations  which  have  the 
same  matrix  T. 

The  method  uses  the  Choleski  decomposition 
T  =  LLr. 

Versions:  MA26A,  MA26AD. 

Language:  Fortran  IV  (standard). 

Date:  May  1975. 

Size:  1.7K;  157  cards. 

Origin:  W.R.Owen,  Univ.  of  Queensland,  Australia. 


MA25A 

To  calculate  the  minim  ax  solution  to  a  system  of 
linear  equations  subject  to  simple  bounds  on  the 

solution  variables.  Calculates  the  vector  x  =  (x . }  that 

J  n 

minimizes 


max  X  a„xi  +  b, 
l </<m  >=l  v  1 

subject  to  \Xj\  <g,  ,/=  1,2, ...,«,  where  n  need  not  be 
equal  to  m. 

More  general  conditions  can  be  applied  by  making  a 
transformation  of  the  variables  and  a  bound  on  the 
accuracy  of  the  maximum  is  returned. 

The  method  is  described  in  Rowell,  Harwell  report 
CSS  11,  (1974)  and  is  based  on  the  exchange 
algorithm  with  a  modification  which  can  short-cut 
some  iterations. 

Versions:  MA25A,  MA25AD. 

Language:  Fortran  IV. 

Date:  December  1974. 

Size:  7.5K;  370  cards. 

Origin:  MJ.D.Fowell,  Harwell,  t 


MA26A 

To  solve  a  system  of  symmetric  positive-definite 
tri-diagonal  linear  equations.  Given  a  matrix 
T  =  Uff}„xn  where  1^=0  when  i>j+l  or/ >/+l,  solves 
the  equations 


X  f nxj  —  bi  l>2,...,n 

J=\  v  J 

The  matrix  is  stored  in  a  compact  form  and  the 
subroutine  may  be  re-entered  with  additional  right 

t  No  longer  at  A.E.R.E.  Harwell. 


MA28A 


To  solve  a  sparse  system  of  linear  equations.  Given  a 
sparse  matrix  A  =  {a,,i  _  this  subroutine 
decomposes  A  into  factors,  solves  Ax  =  b  (or 
optionally  Arx  =  b  ).  It  will  decompose  a  new  matrix 
having  the  same  sparsity  pattern  as  a  previous  one  by 
using  the  same  pivotal  sequence  taking  much  less 
processing  time  than  the  original  factorization.  The 
matrix  A  is  also  allowed  to  be  singular  or  rectangular. 

The  method  is  a  variant  of  Gaussian  elimination  for 
sparse  systems  and  further  information  is  given  in 
Duff,  Harwell  report  R.8730,  (1977). 

Remark:  Supersedes  MAI 8 A. 

Versions:  MA28A,  MA28AD. 

Calls:  MA30A,  MC20A,  MC22A,  MC23A,  MC24A. 
Language:  Fortran  IV  (standard  available). 

Date:  April  1977. 

Size:  9.6K;  653  cards. 

Origin:  I.S.Duff,  A.E.R.E.,  Harwell. 


MA29A 

To  factorize  a  real  symmetric  matrix  and  optionally 
solve  the  corresponding  system  of  linear  equations.  The 
matrix  need  not  be  positive-definite.  Given  a  matrix 
A =  {“tjhxn  there  are  entries  to  factorize  A 

(permuted)  into  the  form  LDLr,  then  given  a  factored 
A,  to  solve  the  system  of  linear  equations 

n 

j^lavxj=bi  1=1,2 . 

by  forward  and  back  substitution,  and  also  an  entry  to 
perform  both  of  these  steps  together. 

The  matrix  is  stored  in  a  compact  form  taking  full 
account  of  symmetry. 

The  method  is  described  in  Fletcher,  'Factorizing 
Symmetric  Indefinite  Matrices',  J. Linear  Algebra  & 
Applies.,  14,  3,  (1976). 


Versions:  MA29A,  MA29AD. 
Language:  Fortran  IV. 

Date:  August  197S. 

Size:  4.7K;  306  cards. 

Origin:  R. Fletcher,  Dundee  Univ. 


MA30A 

Given  a  sparse  matrix  A  =  {Oy)nxn  *  performs  the 
LU  decomposition  of  the  diagonal  blocks  of  a  specified 
permutation  of  A  and  solves  the  linear  equations 

Ax  =  b  ,  i.e. 


or  optionally  solves  Arx  =  b . 

After  decomposing  a  matrix  it  is  possible  to 
decompose  other  matrices  with  the  same  sparsity 
pattern  with  significant  savings  in  processing  time. 
Also,  singular  and  rectangular  matrices  can  be  dealt 
with. 

The  method  is  a  variant  of  Gaussian  elimination  for 
sparse  matrices  see.  Duff  and  Reid,  Harwell  report 
CSS  48  (1977),  and  DufT,  Harwell  report  R.8730 
(1977). 

Remark:  MA28A  provides  a  simple  user-interface  to 
MA30A. 

Versions:  MA30A,  MA30AD. 

Language:  Fortran  IV  (standard  available). 

Date:  June  1977. 

Size:  16.9K;  1186  cards. 

Origin:  I.S.Duff,  A.E.R.E.,  Harwell. 


MB01C 

Given  a  real  matrix  A  =  ia.J  evaluates  the 

(/  nXH 

inverse  matrix  A  1 . 

The  method  is  simple  Gaussian  elimination  with  row 
interchanges.  No  checks  on  the  condition  of  the  matrix 
are  made  and  no  indication  of  the  accuracy  of  the 
inverse  is  given.  If  the  user  is  concerned  about  such 
matters  he  should  consider  the  alternative  subroutines 
in  the  library. 

Remark:  supersedes  MB01B  (removes  the  restriction 
on  the  order  n). 

Version*:  MB01C,  MB01CD. 

Calls:  MC03AS. 

Language:  Fortran  IV. 


Date:  February  1973. 

Size:  2.6K;  91  cards. 

Origin:  MJ.Hopper,  A.E.R.E.,  Harwell. 


MC13D 

Given  the  pattern  of  nonzeros  of  a  spans  matrix  A , 
finds  a  symmetric  permutation  that  makes  the  matrix 
block  upper  triangular,  i.e.  finds  P  such  that 
U  =  PAP  ' 1  is  block  upper  triangular. 

The  method  is  that  of  Tarjan,  SIAM  J.  Computing, 
Vol  1,  (1972),  and  is  also  described  by  Duff  and  Reid, 
Harwell  report  CSS  29  (1976). 

Venions:  MC13D. 

Language:  Fortran  IV  (standard  available). 

Date:  March  1976. 

Size:  2.2K;  131  cards. 

Origin:  I.S.Duff,  A.E.R.E.,  Harwell. 


MC16A 

Appends  an  n+ 1  vector  to  an  nxn  triangular  matrix 

to  form  an  (n+ 1  )x(n+ 1)  triangular  matrix,  i.e.  given  an 
upper  triangular  matrix  U  —  Inx„  and  a  vector 

H  =  (jurM2 . ^n+1 }  this  subroutine  forms  the 

triangular  matrix 


Both  0  and  U  are  stored  in  compact  form. 

Remark:  the  matrix  storage  format  is  that  of  MCI  1  A. 

Versions:  MCI 6A,  MC16AD. 

Language:  Fortran  IV. 

Date:  September  1974. 

Size:  .7K;  29  cards. 

Origin:  MJ.D.FoweH,  A.E.R.E.,  Harwell,  t 


MC17A 

To  delete  a  column  from  an  nxn  triangular  matrix  to 
get  V  =  and  return  a  related  triangular 

matrix  U  =  i>x(n— i)  such  =  VTV. 

Both  the  original  matrix  and  U  are  stored  in  a 
compact  form. 


t  No  longer  at  A.E.R.E.  Harwell. 
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Remark:  the  matrix  storage  format  is  that  of  MC 1 1  A. 

Versions:  MCI 7 A,  MCI 7 AD. 

Calls:  MC11A. 

Language:  Fortran  IV. 

Date:  September  1974. 

Size:  .9K;  38  cards. 

Origin:  MJ.D.Powell,  A.E.R.E.,  Harwell.! 


MC18A 

Given  a  real  symmetric  matrix  A  computes  the 
Householder  transformation  QAQ T  =  T  where  T  is 
symmetric  tri-diagonal.  Also  given  a  vector  x  it  will 
form  Qx  or  optionally  Q  Tx  . 

The  method  is  described  in  Ralston  and  Wilf, 
‘Mathematical  Methods  for  Digital  Computing’,  Vol.  1. 

Versions:  MCI 8 A,  MCI 8 AD. 

Calls:  FM02AS. 

Language:  Fortran  IV  (standard). 

Date:  May  1973. 

Size:  IK;  112  cards. 

Origin:  W.R.Owen,  Univ.  of  Queensland. 


Versions:  MCI 9A,  MCI 9 AD. 

Language:  Fortran  IV  (standard  available). 

Date:  March  1977. 

Size:  2.8K;  167  cards. 

Origin:  J.K.Reid,  A.E.R.E.,  Harwell. 

MC20A 

Sorts  the  nonzeros  of  a  sparse  matrix  from  arbitrary 
order  to  an  ordering  of  columns  or  complete  ordering 
of  columns  and  rows. 

The  subroutine  allows  users  with  sparse  problems  to 
input  the  nonzeros  in  a  natural  order  more  suited  to  the 
individual  problem.  The  resulting  ordered  matrix  can 
subsequently  be  presented  to  other  sparse  matrix 
subroutines  in  the  library. 

An  in-place  sort  algorithm  is  used  which  handles 
each  item  to  be  sorted  exactly  three  times,  so  the 
number  of  operations  for  the  method  is  of  the  order  of 
the  number  of  nonzeros. 

Versions:  MC20A,  MC20AD. 

Language:  Fortran  IV. 

Date:  November  1973. 

Size:  1.6K;  103  cards. 

Origin:  J.K.Reid,  A.E.R.E.,  Harwell. 


MC19A 

To  calculate  scaling  factors  for  a  sparse  matrix 
A  =  \aiJ)nxn  ■  They  may  be  used,  for  instance,  to  scale 
the  matrix  prior  to  solving  a  corresponding  set  of  linear 
equations,  and  are  chosen  so  that  the  scaled  matrix  has 
its  non-zeros  near  to  unity  in  the  sense  that  the  sum  of 
the  squares  of  the  logarithms  of  the  non-zeros  is 
minimized.  The  natural  logarithms  of  the  scaling 
factors  (r(, Cy  ij=  1,2,..., n)  for  the  rows  and  columns 
are  returned  so  that  the  scaled  matrix  has  elements 

bU  =  alJtx^ri  +  CP‘ 

The  subroutine  expects  a  square  nxn  matrix  but 
tolerates  rows  or  columns  consisting  entirely  of  zeros. 
Therefore  a  rectangular  matrix  can  be  handled  by 
embedding  it  in  a  square  matrix  whose  additional  rows 
or  columns  consist  entirely  of  zeros. 

The  method  is  described  in  Curtis  and  Reid,  'On  the 
Automatic  Scaling  of  Matrices  for  Gaussian 
Elimination’.  J.  Inst  Maths.  Applies.  (1972),  10,  pp. 
118-124. 


t  No  longer  at  A.E.R.E.  Harwell. 


MC21A 

Given  the  pattern  of  nonzeros  of  a  sparse  matrix  this 
subroutine  attempts  to  find  a  row  permutation  that 
makes  the  matrix  have  no  zeros  on  its  diagonal. 

The  method  used  is  a  simple  depth  first  search  with  a 
look  ahead  and  is  described  by  Duff,  Harwell  report 
CSS  49  (1977). 

Versions:  MC21A,  MC21AD. 

Language:  Fortran  IV  (standard  available). 

Date:  April  1977. 

Size:  2.6K;  139  cards. 

Origin:  I.S.Duff,  A.E.R.E.,  Harwell. 


MC22A 

Given  a  spans  matrix  A  =  ) _ and  a  row 

\j  nxn 

permutation  matrix  P  and  a  column  permutation 
matrix  Q,  this  subroutine  performs  the  permutation 
A  = PAQ . 

The  nonzero  elements  of  A  are  stored  by  rows  in  a 
compact  form  and  the  user  defines  the  permutation 
matrices  P  and  Q  by  index  vectors  of  length  n. 


i 
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Versions:  MC22A,  MC22AD. 

Language:  Fortran  IV  (standard  available). 
Date:  November  1976. 

Size:  1.3K;  78  cards. 

Origin:  I.S.DufF,  A.E.R.E.,  Harwell. 


MC23A 

Permutes  a  sparse  matrix  to  a  block  lower  triangular 

firm,  i.e.  given  a  sparse  matrix  A  =  {a y)„x„  it 
attempts  to  find  permutation  matrices  P  and  Q  such 
that  X  =  PAQ  is  block  lower  triangular  and  it  returns 
X  to  the  caller. 

Versions:  MC23A,  MC23AD. 

CaOs:  MCI  3D,  MC21A. 

Language:  Fortran  IV  (standard  available). 

Date:  April  1977. 

Size:  3.2K;  202  cards. 

Origin:  I.S.Duff,  A.E.R.E.,  Harwell. 


MC24A 

To  obtain  an  estimate  of  the  largest  element 
encountered  during  Gaussian  elimination  on  a  sparse 
matrix  A  =  i aIJ'lnx„  • The  subroutine  is  given  the  LU 
decomposition  factors  of  A  obtained  from  the 
elimination.  If  the  matrix  has  been  scaled  this  estimate 
will  give  an  indication  of  the  numerical  accuracy  of  the 
decomposition. 

The  estimate  is  obtained  by  using  Hadamard’s 
inequality  on  the  expression  for  the  matrix  elements  in 
terms  of  the  LU  factors  of  the  decomposition.  Further 
details  are  given  in  Erisman  and  Reid,  ‘Monitoring  the 
Stability  of  the  Triangular  Factorization  of  a  Sparse 
Matrix’,  Numer.  Math.,  22  (1974). 

Versions:  MC24A,  MC24AD. 

Language:  Fortran  IV  (standard  available). 

Date:  March  1977. 

Size:  .9K;  46  cards. 

Origin:  I.S.Duff,  A.E.R.E.,  Harwell. 


NB02A 

To  Ibid  a  real  zero  of  a  continuous  real  function  y(x) 

in  a  given  interval  ,  i.e.  solve  the  nonlinear 

equation  in  one  variable 

*0  =  0 


The  limits  a  and  b  must  be  chosen  so  that 
stgn{y(a)}  #  sign{y(b)}  and  the  zero  is  then  located  by  a 
combination  of  linear  interpolation  and  binary 
subdivision.  The  user  must  supply  code  to  evaluate**) 
at  any  point  in  the  interval  and  a  limit  can  be  specified 
for  the  number  of  function  evaluations.  The  subroutine 
attempts  to  refine  the  bracket  down  to  two  points 
where  the  accuracy  parameter  e  is  specified 
by  the  user. 

Remark:  similar  to  NB01A  except  the  convergence 
criterion  is  based  on  x  not  *x). 

Versions:  NB02A,  NB02AD. 

Language:  Fortran  IV. 

Date:  December  1975. 

Size:  1.2K;  113  cards. 

Origin:  A.R.Curtis,  A.E.R.E.,  Harwell. 


OB02A 

Provides  a  simple  graph  drawing  facility.  It  is 
intended  for  very  simple  graphs  and  makes  GHOST 
easier  to  use  by  standardizing  options  so  that  the  user 
has  less  to  remember,  while  still  offering  some  choice. 

There  are  four  options 

(a)  take  a  new  page  and  draw  a  pair  of  labelled  axes 
(restricted  to  the  positive  quadrant). 

(b)  plot  one  of  eight  different  symbols  at  a  specified 
point. 

(c)  connect  a  sequence  of  points  by  straight  lines,  a 
kind  of  curve  drawing  option. 

(d)  print  titles  below  the  graph. 

The  subroutine  uses  the  GHOST  graphics  package 
but  mixing  of  GHOST  subroutine  calls  and  OB02A  is 
not  recommended. 

Versions:  OB02A. 

Calls:  ICO  IAS,  and  various  GHOST  subroutines. 
Language:  Fortran  IV. 

Date:  November  1974. 

Size:  3.8K;  159  cards. 

Origin:  A.R.Curtis  and  MJ.Hopper,  A.E.R.E., 
Harwell. 

OB13A 

To  plot  the  part  of  a  conic  section  that  is  inside  a 
given  triangle.  The  conic  section  is  defined  by  the 
equation  Q(xj>)  =  rj,  where  rj  is  the  contour  height  and 
where  Q(xy)  is  a  quadratic  function  that  is  specified  by 


£ 


its  value  at  the  vertices  and  midpoints  of  the  triangle 
sides. 

The  subroutine  uses  the  GHOST  graphics  package 
but  communicates  with  GHOST  through  only  one 
subroutine  which  makes  it  easily  adapted  for  other 
graphics  packages.  The  method  is  described  in  Marlow 
and  Powell,  ‘A  Fortran  Subroutine  for  Plotting  the  Part 
of  a  Conic  that  is  Inside  a  Given  Triangle',  Harwell 
report  R.8336  (1976). 

Versions:  OBI 3 A. 

Calls:  SHVECS  (a  GHOST  subroutine). 

Language:  Fortran  IV. 

Date:  March  1976. 

Size:  8-4K;  370  cards. 

Origin:  S.Marlow,  A.E.R.E.,  Harwell. 


OBI4A 

To  draw  contours  of  a  Ruction  fixy)  whose  values 
are  specified  over  an  nxm  regular  rectangular  grid. 
Specifically,  the  values  of fixy)  must  be  specified  at  the 
grid  points  (x,  +iAx,  y  +jAy)  0</<m-l,  0</<n-l, 
and  where  Ax  and  Ay  define  the  uniform  spacing  of  the 
grid.  The  contours  at  the  boundary  of  the  rectangular 
region  can  be  improved  if  the  user  can  provide  extra 
values  of  J[xy)  immediately  outside  the  region,  i.e. 
along  any  or  all  of  the  lines  x=xs~ Ax,  x=xs+mJx, 
y=ys~Ay  and y-ys  +nAy. 

The  subroutine  uses  the  GHOST  graphics  package 
but  draws  only  the  contours  and  does  not  draw  axes  or 
perform  mapping  functions.  It  communicates  with 
GHOST  through  the  library  subroutine  OB  13  A  and  is 
easily  adapted  for  other  graphics  packages. 

The  method  is  described  in  Powell,  Harwell  report 
TP  531,  (1973),  and  Marlow  and  Pawell,  Harwell 
report  R.8336,  (1977). 

Versions:  OB  14A. 

Calls:  OB  13 A. 

Language:  Fortran  IV. 

Date:  August  1977. 

Size:  8.3K;  557  cards. 

Origin:  S.Marlow  and  MJ.D.POweUf,  A.E.R.E., 
Harwell. 


OE02A 

To  insert  extra  statements  in  a  Fortran  source  deck 

to  eount  the  number  of  times  state  meats  are  executed 


+  No  longer  at  A.E.R.E.  Harwell. 


during  a  run.  The  user  can  insert  controls  into  the 
source  to  switch  the  counting  on  and  off  at  selected 
points  and  also  request  a  print  out  of  the  counts.  The 
Fortran  is  assumed  to  be  either  standard  ot  IBM 
Fortran  IV. 

To  use  OE02A  an  additional  job  step  is  executed  in 
which  OE02A  will  read  in  the  user’s  source,  insert  the 
extra  statements  and  then  pass  the  modified  source  on 
to  the  next  job  step  which  will  usually  be  the  compile 
step  of  a  compile-and-go  procedure. 

Versions:  OE02A. 

Language:  Fortran  IV  (standard  available). 

Date:  March  1976. 

Size:  14.9K;  426  cards. 

Origin:  J.K.Reid,  A.E.R.E.,  Harwell. 


OE04A 

Allows  a  single  Fortran  deck  to  contain  several 
versions  of  a  program  (or  subroutine)  by  holding 
alternative  statements  as  special  comment  cards,  e.g. 
single  and  double  precision  versions  of  a  subroutine 
which  would  differ  only  in  statements  involving 
precision  specifications. 

OE04A  is  in  essence  a  preprocessor  which  can 
perform  either  a  reversible,  or  a  non-reversible 
transformation  of  one  version  of  a  program  to  another. 

Versions:  OE04A. 

Language:  Fortran  IV  (standard  available). 

Date:  January  1976. 

Size:  9.7K;  138  cards. 

Origin:  J.K.Reid,  A.E.R.E.,  Harwell. 

PE09A 

Given  a  polynomial 

P(x)  =  X  «,<?,(*) 
j=o  J  J 

expressed  as  a  sum  of  orthogonal  polynomials  q/x) 
over  a  point  set  x(/=l,2 . n  this  subroutine  computes 

values  of  P{x),  and  —  at  any  point  x. 

**  dxl 

The  method  is  that  given  by  Smith,  Math.  Comp., 
19,  (1965). 

Remark:  similar  to  PE07A  but  provides  derivative 
values. 
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TB06A 


Vcnlom:  PE09A,  PE09AD. 

Language:  Fortran  IV  (standard). 

Date:  May  1975. 

Size:  .9K;  97  cards. 

Origin:  W.R.Owen,  Univ.  of  Queensland. 


QG03A 

Given  a  piecewise  linear  function  X(x),  defined  by 
points  and  function  values 

Xi=X(il),  /=l,2,...,n,  n>2 

and  given  integration  limits  a  and  b,  and  parameters  tfi 
and  x  the  subroutine  evaluates 

Q(<Pri)=  Ce-^mdx 
Ja 

The  subroutine  can  be  applied  to  the  problem  of 

folding  a  Gaussian  into  a  function. 

Versions:  QG03A,  QG03AD. 

Language:  Fortran  IV. 

Date:  May  1975. 

Size:  2K;  140  cards. 

Origin:  S.Marlow,  A.E.R.E.,  Harwell. 


QG04A 

Given  a  cubic  spline  s(x)  defined  by  knots  4r 
function  values  sf  =  s(4/)  and  derivative  values 

gt  =  jf'sOv)'  1=1.2 2  and  given  integration 

limits  a  and  b,  and  parameters  <p  and  x  the  routine 
evaluates 

£>(<M)=  fde~l<pl(x~xl2s(x)dx 

The  subroutine  is  not  restricted  to  cubic  splines  and 
may  be  used  with  any  piecewise  cubic  function  that  has 
a  continuous  first  derivative  and  which  can  be  specified 
by  its  values  and  first  derivative  values  at  the  joins. 

The  subroutine  can  be  applied  to  the  problem  of 
folding  a  Gaussian  into  a  function. 

Versions:  QG04A,  QG04AD. 

CaOs:  QG02A. 

Language:  Fortran  IV. 

Date:  May  1975. 

Size:  2.1K;  137  cards. 

Origin:  S.Marlow,  A.E.R.E.,  Harwell. 


This  subroutine  computes  tbe  spline,  of  specified 
degree  and  specified  knots,  which  interpolates  a 
arbitrarily  prescribed  function  values.  Specifically, 
given  function  values  /( ,...,/„  at  the  n  data  points 
xt<x2<  ■•■<xfl,  and  given  n-k  knots  i\t, 
/=1,2,...,«-Jt,  (1  <fc<n),  in  the  open  interval 
Xj<x<xn,  this  subroutine  computes  the  coefficients 
a[  of  the  unique  spline  function 

S{x)  =  £  o.N.Xx) 

(=1 

of  degree  k—\  which  has  the  knots  \r).)  and  which 
satisfies  the  interpolation  conditions 
n 

AT.  .(*  )=/.,  j=  1,2 . rt. 

/=  l  ’  J  1 

N.B.  The  function  Nk  j  (x)  is  a  normalized  B-spline  of 
degree  k—  1. 

As  a  by  -product  of  the  calculation  of  the  coefficients 

a] . a  ^  the  subroutine  also  computes  and  returns  the 

value  of  the  integral 

T"  S(x)  dx 

Jx, 

Versions:  TB06A,  TB06AD. 

Language:  Fortran  IV  (standard). 

Date:  December  1976. 

Size:  3.4K;  253  cards. 

Origin:  P.W.Gaffney,  A.E.R.E.,  Harwell.t 


TB07A 

Computes  the  optimal  interpolation  formula  f(x),  of 
specified  degree,  which  interpolates  n  prescribed 
function  values.  Specifically,  given  function  values 
/, ,/2 ,■••,/„  at  the  n  data  points  x,<x2<...<xn  this 
subroutine  computes  the  knots  ,...,9 and  the 
coefficients  al,...,an  of  the  spline  function 

n 

S(x)=  la,NkAx) 

/=! 

of  degree  it— 1,  which  satisfies  the  interpolation 
conditions 

n 

Hatbfkj(Xj)=fj,  j=  1,2 . n 


+  No  longer  at  A.E.R.E.  Harwell. 
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and  which,  whenever  l/^l,  1  is  bounded,  and 
the  value  of  the  bound  is  unknown,  provides  the 
smallest  possible  value  of  B(x)  in  the  error  bound 

IJW  -  £(x)  <  B(x)  max  l/k)(x)l. 

N.B.  The  function  Wfc/(x)  is  a  normalised  B-spline  of 
degree  k—  1. 

As  a  by-product  of  the  calculation  the  subroutine 
also  computes  the  value  of  the  integral 

fXnS(x)dx. 

xi 

The  method  is  described  in  the  Harwell  reports 
Gaffney,  CSS  52  (1977)  and  Gaffney,  R.8781  (1977). 

Versions:  TB07A,  TB07AD. 

Calls:  TB06A,  TB08A. 

Language:  Fortran  IV  (standard). 

Date:  June  1977. 

Size:  1.8K;  103  cards. 

Origin:  P.W.Gaffney,  A.E.R.E.,  Harwell.f 


TB08A 

To  compute  the  knots  of  the  optimal  spline 
interpolation  formula.  Specifically,  the  function  S(x) 
which  interploates  values  of  a  function  of  one  variable 
A x)  at  n  distinct  points  x}<x2<...<xn,  and  which, 
whenever  the  fcth  derivative  of  Ax)  is  bounded  and  the 
value  of  the  bound  is  unknown,  provides  the  smallest 
possible  value  of  B(x)  in  the  error  bound 

| Ax)-S(xXB(x)  max  |/*}(x)!. 

*!<*<*„ 

S(x)  is  a  unique  spline  function  of  degree  k—  1  with  n—k 
knots  Tiv...,rin_k  . 

The  method  is  described  in  the  Harwell  reports 
Gaffney,  CSS  52  (1977)  and  Gaffney,  R.8781  (1977). 

Versions:  TB08A,  TB08AD. 

Language:  Fortran  IV  (standard). 

Date:  June  1977. 

Size:  8K;  550  cards. 

Origin:  P.W.Gaffney,  A.E.R.E.,  HarweU.f 


TG03A 

Evaluates  a  spline  S(x)  of  degree  k—l  and  Its 
derivatives  using  the  B-spline  representation  of  S(x). 


Specifically,  given  knots  and  coefficients  a. ,.../! 

]  fff+a 

in  the  representation 

m+k 

S{x)=  X  «>0,  1, 

(=1 

this  subroutine  computes  the  values  of 
dJ~x 

— — S(x),  j=l,2,...,r,  r^k 
dxJ 

at  a  specified  point  x. 

The  method  is  based  on  De  Boor,  ‘On  Calculating 
with  B-spUnes’,  J.  App.  Theory,  6,  (1972). 

Versions:  TG03A,  TG03AD. 

Language:  Fortran  IV  (standard). 

Date:  January  1977. 

Size:  2K;  131  cards. 

Origin:  P.W.Gaffney,  A.E.R.E.,  HarwelLf 


TG04A 

Given  the  n  points  x,  <x2  <. .  ,<xn  and  a  value  of  x, 
this  subroutine  computes  the  values  of  the  k  B-spUnes 
of  degree  Ar-1 

Mk/*)  1  <Kn-k, 

which  have  knots  at  the  points  xrx/+}  ,...,xj+k  and 
which  cover  the  point  x,  and  it  also  computes  the 
values  of  the  corresponding  integrals 

f*Mkt(Qdt,  KKn-k. 

The  subroutine  uses  the  recurrence  relation  due  to 
De  Boor,  ‘On  Calculating  with  B-splines’,  J.  Approx. 
Theory,  6,  (1972),  to  evaluate  the  B-spline.  The 
integrals  are  evaluated  using  the  method  due  to 
Gaffney,  Harwell  report  CSS  10,  (1974). 

Versions:  TG04A,  TG04AD. 

Language:  Fortran  IV  (standard). 

Date:  December  1976. 

Size:  3.3K;  235  cards. 

Origin:  P.W.Gaffney,  A.E.R.E.,  HarwelLf 


VA13A 

To  calculate  the  minimum  of  a  general  function  of 
several  variables  when  values  of  the  derivatives  with 


f  No  longer  at  A.E.R.E.  Harwell. 
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respect  to  the  variables  can  be  provided,  that  is,  find 
x  =  Uj ,  x2 ,. . .,  xn }  to  minimize  a  function  P( x)  given 


The  subroutine  uses  the  BFGS  variable  -  metric 
method  without  line  searches  of  the  type  analysed  by 
fttwell,  Harwell  report  CSS  15,  (1975). 

Versions:  VA13A,  VA13AD. 

Calls:  MC11A. 

Language:  Fortran  IV. 

Date:  May  1975. 

Slze:4.1K;  180  cards. 

Origin:  M.J.D. Powell,  A.E.R.E.,  Harwell. t 


VA14A 

To  calculate  the  minimum  of  a  general  fanction  of 
many  variables  when  values  of  the  derivatives  with 
respect  to  the  variables  can  be  provided,  that  is,  find 
x  =  I*! , x2,..., xn }  to  minimize  a  function  F(x)  given 


The  subroutine  is  intended  for  problems  that  are  so 
large  that  an  nxn  matrix  cannot  be  stored 
conveniently. 

A  conjugate  gradient  method  is  used  that  includes 
an  automatic  restart  procedure  and  it  is  described  in 
POwell,  Harwell  report  CSS  24,  (1975). 

Versions:  VA14A,  VA14AD. 

Language:  Fortran  IV. 

Date:  November  1975. 

Size:  3.8K;  189  cards. 

Origin:  MJ.D.Powell,  A.E.R.E.,  Harwell.f 


VD03A 

Given  a  curve  and  a  data  point  this  subroutine  will 
compute  the  point  on  the  curve  which  Is  nearest,  in  a 
weighted  least  squares  sense,  to  the  data  point  That  is, 
given  a  function  y(x)  representing  the  curve  and  a  point 
(ij?)  with  associated  weights  w  and  w  ,  the 
subroutine  minimises 

S(x)  =  xf  +  yf 

Values  of  the  first  and  seeond  derivates  of  y{x)  must  be 
supplied  by  the  user. 


The  method  is  Newton-Raphson  with  a  bound  on 
the  x  value  to  ensure  convergence. 

Versions:  VD03A,  VD03AD. 

Language:  Fortran  IV  (standard). 

Date:  November  1974. 

Size:  2.6K;  342  cards. 

Origin:  W.R.Owen,  Univ.  of  Queensland. 


ZA19AS 

Allows  a  Fortran  program  to  disable  or  enable 
exponent  underflow  interrupts.  The  old  status  is  always 
returned  to  the  user,  thus  allowing  him  to  reset  to  a 
previous  condition  if  desired. 

Fortran  error  processing  of  exponent  underflows  can 
be  expensive.  If  a  program  is  beyond  its  developement 
stage  and  many  valid  underflows  are  expected,  which 
are  to  be  fixed  to  a  zero  result,  then  ZA19AS  should  be 
used  to  save  c.p.u.  time. 

Versions:  Z A 19 AS. 

Language:  370/BAL. 

Date:  May  1975. 

Size:  .IK;  91  cards. 

Origin:  MJ.Hopper,  A.E.R.E.,  Harwell. 


ZA21AS 

To  allow  the  user  of  the  IBM  Fortran  DEBUG 
facility  to  switch  the  INIT  option  on  and  off  at  chosen 
points  in  a  program. 

Use  of  ZA21AS  can  reduce  the  amount  of  printed 
output  produced  by  a  debug  run  and  it  also  makes 
some  reduction  in  the  execution  time. 

ZA21AS  is  based  on  IBM  code  published  in  the 
S.E.A.S.  Newsletter,  4,  (June  1975),  under  the  name 
INIT. 

Versions:  ZA21AS. 

Language:  370/BAL. 

Date:  November  1975. 

Size:  .IK;  40  cards. 

Origin:  MJ.Hopper,  A.E.R.E.,  Harwell. 


t  No  longer  at  A.E.R.E.  Harwell. 
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